/* 
 * File:   vecs.h
 * Author: boyan
 *
 * Created on December 12, 2007, 12:56 PM
 */

#ifndef _VECS_H
#define	_VECS_H

#include "glodef.h"
// #include <petscvec.h>

typedef struct _vector * vector;

typedef struct _scalar * scalar;

#ifdef	__cplusplus
extern "C" {
#endif

    scalar vec2scal(vector Vec);
    vector scal2vec(scalar Scal);

    vector alloc_vec(int N);
    void free_vec(vector A);
    void copy_vec(const vector S, vector T); /* copies contents of S into T */
    void clean_vec(vector A); /* sets contents of A to 0 */
    void print_vec(vector A, char *fname);

    void scale_vec(vector y, double a); /* y[i] = a * y[i] */
    void axpy_vec(vector y, double a, vector x); /* y[i] = a * x[i] + y[i]*/
    // void axpby_vec(vector y, double a, double b, vector x); /* y[i] = a * x[i] + b * y[i]*/
    void waxpy_vec(vector w, double a, vector x, vector y); /* w[i] = a * x[i] + y[i] */
    double dot_vec(vector x, vector y); /* return x'. y */

    // void get_ownership_range_vec(vector A, int *low, int *high);
    void set_value_vec(vector A, int pos, double *val);

    int get_N_vec(vector A);
    void check_use_vector(vector A);
    double *get_XYZ_vec(vector A, double **X, double **Y, double **Z);
    void restore_XYZ_vec(vector A, double **X, double **Y, double **Z);
#define use_vector(A)     {  double *A ## X=0, *A ## Y=0, *A ## Z=0; \
            if(A) { check_use_vector(A); \
            get_XYZ_vec(A, &A ##X, &A ## Y, &A ## Z); }
#define done_vector(A) if(A) restore_XYZ_vec(A, &A ##X, &A ## Y, &A ## Z); }

    scalar alloc_scal(int N);
    void free_scal(scalar A);
    void copy_scal(const scalar S, scalar T); /* copies contents of S into T */
    void clean_scal(scalar A); /* sets contents of A to 0 */
    void print_scal(scalar A, char *fname);

    void scale_scal(scalar y, double a); /* y[i] = a * y[i] */
    void axpy_scal(scalar y, double a, scalar x); /* y[i] = y[i] + a * x[i] */
    void waxpy_scal(scalar w, double a, scalar x, scalar y); /* w[i] = a * x[i] + y[i] */
    double dot_scal(scalar x, scalar b);
    double sum_scal(scalar A);

    void get_ownership_range_scal(scalar A, int *low, int *high);
    void set_value_scal(scalar A, int pos, double val);
    void assemble_scal(scalar A);

    int get_N_scal(scalar A);
    void check_use_scalar(scalar A);
    // void *get_Vec_scal(scalar A);
    double *get_F_scal(scalar A, double **F);
    void restore_F_scal(scalar A, double **F);
#define use_scalar(A)    { double *A ## F=0; \
            if(A) { check_use_scalar(A); \
            get_F_scal(A, &A ##F); }
#define done_scalar(A) if(A) restore_F_scal(A, &A ##F); }


#ifdef	__cplusplus
}
#endif

#endif	/* _VECS_H */

